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I. INTRODUCTION 

With gravitational wave detectors such as LIGO, 
VIRGO and GEO operational the problem of calculat- 
ing gravitational wave templates has become even more 
urgent. Amongst the physical models that are the best 
candidates for producing detectable wave signals those 
including highly relativistic matter near black holes stand 
out. In cases such as black hole / neutron star binaries, 
binary neutron star systems that collapse promptly to 
a black hole, accretion flows onto black holes and many 
models for gamma-ray bursts, detailed numerical simula- 
tions will be required to find the impact of varying phys- 
ical parameters on the gravitational waves produced. 

Black hole / neutron star binaries are on astrophys- 
ical grounds believed to be as likely as binary neutron 
star mergers, with expected event rates of one per year 
in a sphere of about 70 Mpc radius 0. While signals 
from binary neutron stars are expected to give us in- 
formation about the masses, spins and locations of the 
objects, they are not expected to give information about 
the internal structure of the stars. Signals from mixed 
binary systems, on the other hand, will provide infor- 
mation about the neutron star structure and equation of 
state (EOS) 0. 

The crucial problem for numerical simulations involv- 
ing 3D general relativity which must be overcome to sim- 
ulate such physical systems is stability. With current 
formulations of the vacuum Einstein equations it is pos- 
sible to produce long term simulations of black holes in 
certain situations These simulations typically 

require some part of the computational domain inside the 
black hole to be excised, with an inner boundary condi- 
tion placed on a surface inside the apparent horizon. This 
apparent horizon is never outside the event horizon. Be- 
cause no physical signal can travel outwards from such an 
horizon, excising the interior (or parts of it) should not 
affect the exterior spacetime, which is the only region 



we can observe. The main reason to excise the part of 
the spacetime containing the singularity is that otherwise 
steep gradients near the physical singularity form, which 
numerical codes cannot handle. The only long term sta- 
ble simulations including black holes performed without 
excision 0, Q use techniques that are only applicable 
when the black hole is present in the initial slice. 

In contrast most simulations including hydrodynamics 
have either been performed on a fixed spacetime back- 
ground, or have only been run until a short time after 
the formation of the black hole There have been 
fully dynamical simulations of matter with black holes in 
axisymmetry such as , but few in 3D 0, . 

In this paper we will present a simple method for 
excision boundaries applied to hydrodynamics. The 
boundary condition is based on High-Resolution Shock- 
Capturing (HRSC) methods which may be used in a hy- 
drodynamics code, and theoretically could be applied to 
any system using such HRSC methods. We show how 
it can be applied to three different, standard reconstruc- 
tion schemes: TVD, ENO and PPM. Because of some 
problems found using this boundary condition with the 
PPM scheme we introduce a modified version of PPM 
(MPPM), which solves these problems. This excision 
method, combined with a suitable excision method for 
the spacetime, allows long term simulations of matter in 
black hole spacetimes. This has been shown e.g. using 
our hydrodynamics code called Whisky in lOj. 

The outline of this paper is as follows. In section lU 
and ^Q] respectively we outline the equations and HRSC 
methods that will be used. The modifications required 
at excision boundaries are given in section Hvl Section Ivl 
contains the tests used to validate the boundary con- 
ditions. Throughout this paper we shall use geometric 
units where c = G — Mq = 1. Greek indices are taken 
to run from to 3, latin indices from 1 to 3. We adopt 
the standard convention for the summation over repeated 
indices. In section iHll latin indices denote the cell index. 
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II. MODEL AND EQUATIONS 

We are interested in simulating hydrodynamical flows 
near black holes, so we use the equations of general rel- 
ativistic hydrodynamics coupled to a dynamical space- 
time, described by full general relativity. We use the 
flux-conservative Valencia 3-1-1 formulation of the hydro- 
dynamical equatio ns [T^ [T3 . [l5j | . Although the results of 
Lax and Wendroff [1^ show that if the scheme converges 
then it will converge to one of the (possibly infinitely 
many) weak solutions of the system of equations, the use 
of the correct conservation law form is known to be cru- 
cial as the results of Hou and Lefloch [T^ show that a 
non-conservative scheme will generically converge to the 
wrong weak solution. The stress energy tensor is written 
as 

Tf_,i, = phu^Ui, + pg^^, (2.1) 

where p is the mass density of the fluid, h — 1 + e + 
the specific enthalpy with e the specific internal energy, 
and p the pressure, is the 4-velocity of the fluid and 
g^u the 4-metric of the spacetime. 

The equations of relativistic hydrodynamics can then 
be written in a conserved form 

atq-f9,f«(q)=s(q). (2.2) 

Here we use the Valencia form [TlIT^ as in [TollTlligt. 

The q is a vector of conserved variables. They are not 
strictly conserved because of the source terms s, but in 
flat space the sources s vanish. 

The system requires an EOS to close it. This is nor- 
mally given by specifying the pressure p in the form 
p = p{p,e). In this paper we will consider the standard 
perfect fluid gamma-law EOS 

P={T-l)pe, (2.3) 

and often we will restrict to the polytropic EOS 

p = Kp^, (2.4) 

which is a good approximation to neutron star matter in 
a cold neutron star and in the absence of shocks. 



III. NUMERICAL METHODS 

The numerical methods that we use for the evo- 
lutions of the hydrodynamic variables are all High- 
Resolution Shock-Capturing (HRSC) methods. We 
take the semidiscrete or method of lines reconstruction- 
evolution viewpoint: a piecewise continuous interpolant 
(the reconstruction) of each variable is found, at each 
cell boundary a Riemann problem is solved (usually ap- 
proximately), the time derivative for each variable is con- 
structed from the flux differences through the boundaries 
of each cell and the source terms are calculated, and the 



solution at the new timelevel is found using some suitable 
time integrator. 

In what follows we shall specialise to the case of a uni- 
form Cartesian grid. 

A. Riemann solvers 

Our code implements three independent approximate 
Riemann solvers (HLLE, Roe, modifled Marquina). Al- 
though an integral part of the full HRSC method they are 
irrelevant for the problem of excision, as explained below. 
The Riemann solver employed in all tests shown below is 
the modifled Marquina solver described in |20i, |21|| . 

B. Reconstruction methods 

Four separate reconstruction methods are considered 
here. Each are applied in a dimensionally split fashion; 
the reconstructions are performed along each coordinate 
axis in turn. Although some of the reconstruction meth- 
ods considered here have formal orders of accuracy that 
are better than second order, the formal global order of 
accuracy of the code is at best second order. This is due 
to the extension to multiple dimensions and the coupling 
to the spacetime. The flux through the cell boundary is 
approximated by the flux through the point in the mid- 
dle of the cell boundary, using the reconstructed fluid 
variables and a second order approximation of the met- 
ric terms. The use of reconstruction methods that have 
one dimensional formal orders of accuracy better than 
second order leads to improvements in actual accuracy 
but not in the formal convergence order. 

The method with the lowest formal order of accu- 
racy is a slope limited total variation diminishing (TVD) 
method. The essentially non-oscillatory (ENO) methods 
may have extremely high orders of accuracy and are more 
accurate in absolute terms than TVD. However, it is of- 
ten more efficient to use the piecewise-parabolic method 
(PPM) even though it only has at most third-order accu- 
racy. The fourth method is a variant of PPM, and it is in- 
troduced for the first time here. It solves some problems 
found using PPM near the excision region. Therefore we 
caU it MPPM for 'modified PPM'. 

For conventions about subscripts and superscripts to 
denote cells, cell boundaries and left or right sided cell 
boundary values also see Fig. ^ 

1. Slope-limited TVD 

Total variation diminishing (TVD) methods as given 
by, e.g., [i^, ensure that the solution remains monotonic 
to avoid spurious oscillations near discontinuities. In the 
case of slope-limited TVD the reconstruction is given by 
an average of a first-order and a second-order reconstruc- 
tion. To reconstruct the variable q in the cell centred at 
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Xi two local "slopes" are defined, 

A+ = Qi+i - Qt, 



(3.1) 



which are then averaged to get the slope in the cell = 
(A~ + Af)/2. This slope is then multiplied by a limiter 
which is a function cj) of the local slopes (p = ^Pi^T ^ ^t)- 
This limited slope A,; then gives the cell boundary data 
by 



9i±i/2 = g» ± ^A, 



(3.2) 



If the slopes are not limited (i.e., the limiter (p — 1) 
then the method is second order accurate. However, in 
the presence of steep gradients or discontinuities such a 
reconstruction will be oscillatory due to Gibb's effects. 
In these regions the limiter reduces the slopes to avoid 
overshoots, retaining monotonicity of the reconstruction. 



2. ENO 

The ENO methods of Harten et al. [2^ are a very 
general class of methods. Here we only consider the 
simple ENO reconstruction of the variables as given by 
Shu in |24| . This provides arbitrary order of accuracy in 
space. 

As in the TVD case, ENO reconstructions are cell 
based. The p"^ order reconstruction compares all pos- 
sible polynomial stencils of size p containing the cell to 
be reconstructed. The "least oscillatory" stencil is chosen 
by minimizing the absolute values of the divided differ- 
ences that make up the stencil. 



Let p be the order of the reconstruction. Suppose we 
are reconstructing the scalar function q in cell i. We start 
with the cell i. We then add to the stencil cell j, where 
J — J ± 1, where we choose j to minimize the Newton 
undivided differences 



q[i-l,i] = qi- qi-i, 
q[i,i+l] = gj+i - qi. 



(3.3a) 
(3.3b) 



We then recursively add more cells, minimiz- 
ing the higher order Newton divided differences 
q[i — 1, . . . ,i + j] defined by 



q[i-l,...,i+j] = q[i,...,i+j]- 

q[i-l,...,i+j -I] 



(3.4) 



The reconstruction at the cell boundary is given by a 
standard p"' order polynomial interpolation on the cho- 
sen stencil. 

Shu |2^] has outlined an elegant way of calculating the 
cell boundary values solely in terms of the stencil and the 
known data. If the stencil is given by 

S{i) = {i-r,...A+p-r-l}, (3.5) 
for some integer r, then there exist constants c^j depend- 
ing only on the grid Xi such that the boundary values for 
cell li are given by 



qi+i/2 



p-i 



1=0 
~1 

C. 



Ep— 1 



(3.6) 



The constants o, 
formula 



are given by the rather complicated 



This calculation simplifies considerably if the grid is 
evenly spaced. For this case the coefficients up to 
seventh-order are given by Shu [2^. 

3. PPM 

The Piecewise Parabolic Method (PPM) of Colella 
and Woodward |2 a. g eneralized to relativistic flows by 
Martf and Mtiller |2y|, is third-order accurate in space 
with particular special cases to ensure monotonicity at 



shocks, sharpening of contact discontinuities, and shock 
detection. 

The outline of the general PPM method is as follows. 
The first step is to interpolate a quartic polynomial to 
the cell boundary, 

gi+i/2 = ^ {qi+i + *) + ^ (^mft - ^mqi+i) , (3.8) 

where 
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C mm{\qi+i - qi^i\,2\qi+i - qi\,2\qi - qi-i\) sign(gi+i - g^-i) if {qi+i - qi){qi - qi-i) > 
6„,q, = I . (3.9) 

I otherwise 



At this point we set both left and right states at the 
interface to be equal to the interpolated value, 

qt+i/2 = 17+1/2 2- (3-10) 

This reconstruction will be oscillatory near shocks. Be- 
fore a step which will preserve the monotonicity there are 



two other steps, which may be applied. 

Firstly we may "steepen" discontinuities. This is to 
produce sharper profiles and is only applied to discon- 
tinuities that are mostly a contact (see for the de- 
tails). This procedure replaces the cell boundary recon- 
structions of the density with 



Pt-i/2 = Pti/2(1 (^Pt-i + ^S,nPi-i^ V, (3.11a) 



(3.11b) 



P^+l/2 = /'^+l/2(l ' ^) + ^^+1 " V, 

where rj is defined as 

rj = max [0,min (1,771(77 - 772))] , (3.12) 

where 771 , 772 are constants and 

--^ if - d^p,+id^pi^i > and \pi+i - pi^i \ - es mm[\p,+i\, \pt-i\j > 

D Pi+i - Pi-l 

otherwise 

(3.13) 

with es another constant and 

S'^Pi = pi+i - 2pi + pi-i. (3.14) 
I 

Suggested values for the constants 771,772 and eg can be tonicity enforcement is the "flattening" of the zone struc- 
found in ture near shocks. This adds simple dissipation, altering 

Another step that may be performed before mono- the reconstructions to 



g+1/2 = v^q+_^/^ + {l~v,)q,, 

9i+l/2 = V^qi+lJ2 + {'^-vM^^ 



(3.15) 
(3.16) 



where 



max 



0, 1 — max (J, [x'2 ^1 

Pi+2 - Pi~2 



if e min(pi_i,pi+i) < pi+i - pi-i and ?;f_i - 7;f+i > 
otherwise 

I 



(3.17) 



and 0^1,^2 and e are again constants with suggested val- ues given in |25|. Note that this step is a simplification 
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of the one in 25, 26], which reduces the communication The final step is apphed to preserve monotonicity, 
overhead in parallel runs. No problems were encountered which prevents oscillatory reconstruction near shocks, 
using this modification. The following replacements are made: 



4-1/2 = 3<Z, - 2g- if ((77^^/2 - g+ - 1(4-1/2 + 17+1/2)^ > ^(*r+i/2 - 9^^-1/2)"' 



i+1/2 



3q,-2q+2 if {q^^^ 



+1/2 '^i-1/2' 



* - 1(4-1/2 + 1^+l/2)) < -li1^+l 



/2 I1-I/2) 



r 



(3.18a) 
(3.18b) 

(3.18c) 



4. MPPM 

As we shall see below there are situations where the 
PPM scheme does not give good results. This is due to 
the use of a fixed stencil in the first part of the PPM algo- 
rithm where the interpolated value is computed. This can 
cause problems with supersonic flows and when, as in our 
excision method, we try to match to a non-PPM method. 
To avoid this problem the modified PPM scheme changes 
the interpolated value in (|3.8() . The modified method 
uses the same steps as the PPM scheme after the inter- 
face values are set in Eq. H3.10|l (discontinuity steepening, 
zone flattening and monotonicity preservation). 

PPM uses the points at« — l,z,z-t-l,?-t-2for the fourth 
order interpolation function to obtain a value at i -|- i . In 
the case of supersonic flow this centred stencil takes too 
much information from one side and too little from the 
other. Because of this, the method may produce small 
amplitude oscillations that are stable (due to the other 
steps in the PPM method) but which do not converge 
away. To cure this we need to use a stencil that depends 
on the data in some way. 

In particular, we calculate the minimum and max- 
imum characteristic velocities from the minimum and 
maximum eigenvalues A_ and A-|_ of the Jacobian ma- 
trix df/dq. Depending on these values, we define 



A_ +A+ 
lA-l + IA, 



(3.19) 



which is a continuous function of the characteristic speeds 
that has value ±1 whenever the flow is supersonic. Be- 
cause the eigenvalues can only be in the range [—1,1] a 
is also in that range. 

Using a we calculate the initial boundary states 



='4+1/2 



H 9l + (1 - H) qi+1/2, if a < 0, 
l"| gfl + (f - |a|) 91+1/2, if a > 



(3.20) 



with 



<1L 



12 
1 

12 



(13^4+1 - 5q.,+2 + qi+3 + 3gi) , (3.21) 
(13g, - 5g,_i + g,_2 + 3g,+i) . (3.22) 



In this way we shift the stencil at most one cell to the side 
depending smoothly on the data. Note that this choice is 
arbitrary and could be chosen in a different way, but was 
sufficient for the tests used here. Because the value of the 
interpolating function at the boundary can be outside the 
range given by the values at the neighbouring cells qi and 
qi+i, we set the boundary values q^^i/2 ^^'^ 1i+i/2 
case to be equal to the closest nearby cell value. Once the 
interpolated value is set we proceed as in the PPM case 
with the steps required to, e.g., preserve monotonicity, 
as given in Eqs. H3.11H3.18|) . 

This procedure takes the characteristic speeds of the 
data into account and should therefore be more suit- 
able to extreme cases where both minimal and maximal 
speeds have the same sign and are not close to zero. How- 
ever, in cases with minimal and maximal characteristic 
speeds of the same size, but different sign, it should be 
as good as PPM. 



IV. EXCISION BOUNDARIES 

Excision boundaries are based on the principle that a 
closed region of spacetime that is causally disconnected 
from the rest of the simulated spacetime can be ignored 
without affecting the results in the exterior spacetime. A 
black hole event horizon is given indeed by the boundary 
of such a causally disconnected region of spacetime. As 
generically a curvature singularity will form within the 
horizon it will be necessary to remove it, and the associ- 
ated steep gradients, from the numerical domain. 

Excision boundaries are usually placed within a 
trapped surface such as an apparent horizon, as it is not 
possible to find the event horizon locally in time. We 
note that for a cubical region that is excised on a Carte- 
sian grid to be a true trapped surface it may have to be 
placed well within the horizon, as pointed out by |27j| . 
However, inflow boundary conditions applied to surfaces 
within the apparent horizon have been shown to work 
well 0. 

The same considerations apply to the hydrodynamical 
variables as to the spacetime field variables. Again there 
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1 X 1 X 

Copy to set up reconstruction. Force a trivial Only fluxes for the boundaries of cells outside 

Riemann problem at the excision boundary. the excised region are required. 



FIG. 1: A schematic view of the excision algorithm. The excised region is indicated by the shaded grey cells and the excision 
boundary by the vertical dotted line. In the left panel we indicate how the reconstruction method is modified near the excision 
boundary. The values required to compute the fluxes for the cell next to the boundary are indicated, along with the method 
used (simple copying) to "extend" the data across the excision boundary. Using this we can compute the fluxes through the 
cell boundaries for all cells outside the excised region. As shown in the right panel, these are the only fluxes required. 



will be a region of spacetime in which the fluid is causally 
disconnected from the exterior. In practise due to the 
coupling with the spacetime field this will coincide with 
the event horizon as for the spacetime field. However, if 
the spacetime field variables have been excised then in 
principle there is no need to excise the hydrodynamical 
field variables at all, as these are not expected to blow 
up and form singularities independently of the spacetime 
field variables. 

In practise unacceptable numerical errors are found if 
the hydrodynamical variables are not excised. In partic- 
ular at typical resolutions it is possible for numerical ef- 
fects in the interior of the black hole horizon to propagate 
out of the horizon through the hydrodynamical variables. 
Although we will exploit the independence in principle 
of the use of excision boundaries for the different field 
variables later, in practise both should be used simulta- 
neously. 

Excision boundary methods can be viewed in a number 
of different ways. Viewed at the level of the differential 
equations, all characteristics of physical quantities are 
going into the excision region. At the level of the Rie- 
mann problem, all possible waves from the solution of the 
Riemann problem must be contained within the excision 
region. 

This immediately suggests a method that could be 
used at the excision boundary. As a HRSC method nat- 
urally changes the stencil locally depending on the data, 
we just need to guarantee that the data used to construct 
the flux at the excision boundary is the physically rele- 
vant data outside the excision region. This also suggests 
that the natural choice to put the excision boundary is 
not a cell centre, but a cell boundary. 

While very simple, this method is not stable in general. 
It is equivalent to filling the first cell inside the excision 
region by linear extrapolation from the exterior (for a 
second-order HRSC method such as slope limited TVD 
reconstruction). This is not guaranteed to reduce the to- 
tal variation of the solution, and so even simple examples 
do fail with this boundary condition. 



On the other hand the simplest outflow boundary con- 
dition at the excision boundary does not have this prob- 
lem. In particular, we may apply zero-order extrapo- 
lation (a simple copy) to all variables at the boundary 
to create data at the first cell inside the excision region. 
This is done for each of the three reconstruction meth- 
ods. If one method requires more cells for the stencil in 
the interior of the excision region, we can then force the 
stencil to only consider the data in the first cell and the 
exterior as above. 

A summary of our method is as follows: 

1. Any point for which all possible reconstruction 
stencils are contained within the exterior spacetime 
(i.e., do not contain an excised point) is evolved in 
the normal fashion. 

2. Any point which is contained within the excised 
region is not evolved at all. 

3. All other points have at least one stencil that con- 
tains an excision point. To update these points we 
must reconstruct the data to both sides of each cell 
boundary. Once this is done the Riemann problem 
can be solved, the boundary flux found and the up- 
date computed. In order to reconstruct the data we 
do the following: 

(a) The values for all variables in the cell next to 
the excision boundary are copied to the first 
cell within the excision boundary. 

(b) A reconstruction of all cells is produced us- 
ing only the data outside the excision region 
and the first cell within. The precise method 
depends on the global reconstruction method 
used and is explained below. In outline, the 
idea is to retain the order of accuracy by re- 
taining the same size of the stencil if possible. 

(c) This method produces reconstructed data on 
both sides of every cell boundary except for the 
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excision boundary, where only the reconstruc- 
tion from the exterior can be produced. At 
this boundary we assume on physical grounds 
that the solution of the Riemann problem 
must be given by the exterior state (as this is 
assumed to be an outflow boundary). There- 
fore at these points we do not solve the Rie- 
mann problem, but calculate the flux from the 
exterior reconstruction (see Fig.^. 

As an aside, we note that this means the method is 
independent of the choice of Riemann solver. However, 
the choice of a trivial Riemann problem at points that 
are not excised but where the reconstructed states are 
identical (e.g., in the atmosphere) can lead to gains in 
computational efficiency. 

To completely describe our method in practice we just 
need to describe how the stencil is altered with the spe- 
cific reconstruction methods we use. In what follows we 
consider a set of cells in one dimension. The coordinate 
is labelled x, the cell centres x^, and the cell boundaries 
^i±i/2 in an obvious notation. The excision boundary 
will be denoted xb+i/2, with the cell xb+i being excised. 
We want to calculate the update term for the cell centred 
at xb by calculating the fluxes at xb±i/2- 

Excision of the spacetime variables is applied using the 
simple lego excision method described in, e.g., 

1. Slope-limited TVD 

In this case it is clear that only the reconstructions at 
xb±i/2 ai'e affected by the excision region. It is also clear 
that setting 

As = i(A-+A+)=0, (4.1) 

ensures that only data outside the excision region is used 
and is consistent with the TVD reconstruction. 



2. ENO 

As described in section llll B 2l the ENO method uses a 
stencil width depending on the desired order of accuracy 
p. Hence the reconstruction in the cells centred between 
Xb and xb-p+2 are affected by the excision region. 

It is clear how to ensure that the stencil chosen by the 
ENO reconstruction does not include any points from 
inside the excision region. The first divided differences 
that include points within the excision region are set to 
extremely large values with oscillating and growing am- 
plitude such as {—iy x 10^'^ where i is the index of the 
cell {i > 0). The higher order divided differences are 
calculated from these, ensuring that the least oscillatory 
stencil does not include points from the excision region, 
but more points away from it. This is__a simple applica- 
tion of the principle outlined by Shu 24] for dealing with 
outflow boundaries with ENO methods. 



3. PPM 

PPM uses a point stencil of five points, but recon- 
structs around a cell boundary instead of in a cell. Hence 
the reconstructions affected by the excision region are 
the reconstruction at left edge of xb+i/2 and both recon- 
structions at Xb-1/2- 

We have not attempted to find a consistent third-order 
reconstruction for these points. Instead we use the iden- 
tical reconstruction as in the TVD case of section llll B II 
As this provides a second-order stable reconstruction at 
an outflow boundary it does not affect the order of accu- 
racy in the rest of the domain. 

4. MPPM 

Since excision should only be applied on a boundary 
with only outgoing characteristics, a should only be — 1 
or 1 depending on the direction of the boundary. Because 
of the symmetry of the cases, without loss of generality 
we assume a — 1 from here on, that is assuming the 
excision boundary to be to the right of the computational 
domain. 

In the case a = 1 we effectively use a stencil with one 
point to the right and three points to the left instead of 
the usual symmetric four point stencil of the interpola- 
tion step of PPM. Because of this, we can use the usual 
MPPM procedure in the whole domain except for only 
one point next to the excision boundary. There we set 
the left and right boundary values to be the cell value as 
in the TVD case. 



V. NUMERICAL TESTS 

In this Section we present some simple numerical 
tests to validate our approach. The results were cal- 
culated using the Whisky code [l^. Where the space- 
time was evolved an implementation of the NOK 
(BSSN) formulation of |30l . Isij was used, with the exci- 
sion method given in Apparent horizons were found 
using the code described in l32l and event horizons us- 
ing the code described in [33 ■ AH codes use the Cactus 
framework [s^. 

A. Shock wave tests 

Special relativistic shock tubes are the simplest possi- 
ble test. For this test the spacetime background is fixed 
to be Minkowski spacetime in standard coordinates. It 
is simple to choose initial data such that the chosen ex- 
cision boundary is an outflow boundary for all time. In 
particular, we are using initial data in which = 10, 
PR = 1,VL^VR^ 0, PL = 40/3 and pr = 2/3 x 10"^ 
and the velocity is normal to the interface. The ideal 
fluid equation of state, Eq. (|2.3|) . is used, with F = 5/3. 
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We have considered three different cases. The first 
and simplest was a shock ahgned with a coordinate axis 
with the excision boundary normal to the shock. The 
second and numerically more complex was a shock along 
the diagonal of the 3D box with the excision boundary 
normal to the shock. The final test consisted of a shock 
not aligned with the grid and with the excision boundary 
given by a sphere. The test requires more care in setting 
the domain evolved as the excision boundary is no longer 
an outflow boundary over the entire spacetime, but only 
within a localized region. Here we set the domain such 
that the intersection of the excision boundary with the 
domain forms a hemisphere. 

The results of the first test for PPM are shown as an ex- 
ample in Fig. 121 All methods give nearly the same results, 
are stable and for all methods the shock passes cleanly 
through the excision boundary, with the results matching 
the analytic solution well. All other tests showed similar 
results, with stable and accurate results independent of 
the relative geometry of the shock and excision boundary. 

An example of what may go wrong with an incorrect 
excision boundary condition is shown in Fig. 13 Here the 
standard excision boundary condition presented in sec- 
tion ^3 is compared to an alternative method where the 
first point inside the excision region is set to very small 
values. The result is a stable evolution which is correct 
except for a few points near the boundary. Here the den- 
sity overshoots except for the point nearest the bound- 
ary which itself drops to atmosphere values. Although 
in this simple test the boundary condition is inaccurate, 
but stable, this seems unlikely to remain the case in more 
complex situations. Other boundary conditions, such as 
first-order extrapolation to the first point within the ex- 
cision region, fare even worse, as they actually fail to 
produce a stable evolution even on this simple test. 



B. Michel solution 

The Michel solution [s^ is a stationary solution of 
spherical accretion onto a Schwarzschild black hole in 
the test fluid approximation. Here we will consider 
Eddington-Finkelstein coordinates so that the slice pene- 
trates the horizon. This solution is often used to validate 
relativistic hydrodynamic codes 36, 37]. 

The solution is derived from the basic equations of 
mass and energy conservation for the fluid and is ob- 
tained numerically using a Newton-Raphson interation 
scheme. 

The first three plots in Fig.^jshow a part of the results 
of a series of runs made to test the Michel solution. In 
these plots the ENO reconstruction method is used, but 
the graphs look very similar for TVD, PPM and MPPM. 
We set the Michel solution to be our initial data, then 
we evolve it and verify its accuracy. In the plot, the 
exact solution is compared with the computed solution 
at time t = 40 (in units of the background black hole 
mass). Using a logarithmic plot we show the results for 



the density, x-velocity and pressure. The two sets of 
data are in very good agreement, as expected from the 
staticity of the solution. It is possible to notice a slight 
discrepancy in the first grid point next to the excision 
region (dotted line). This discrepancy is related to the 
lower accuracy which we have there, visually amplified 
on the plot by the choice of logarithmic scale. This is the 
case for all reconstruction methods. 

The fourth plot in Fig.Qlshows instead the convergence 
factor obtained for all the reconstruction methods on the 
same set of simulations. As evident in the figure, PPM 
does not show convergence, while the other reconstruc- 
tion methods show the expected convergence factor of 
approximately two. We will discuss the possible reason 
for that later. 



C. Neutron Star collapse 

As the next test we consider the collapse of a neutron 
star, in a simulation in which both the spacetime and the 
matter are evolved. This test is considerably more dif- 
ficult than previous tests, requiring both hydrodynamic 
fields and spacetime to be consistently excised, and re- 
quiring a dynamic excision region. 

As in [l9j a spherically symmetric polytrope with 
K = 100, r = 2 and a central density of pc = 8 x 10"^ « 
4.95 X IQ-^^gcm"^ is used. The collapse was induced by 
lowering K by 2% and rescaling the matter variables to 
enforce the solution of the constraints. Note that be- 
cause this satisfies the constraints and the spacetime is 
not changed, the ADM mass of the system is also not 
changed by this. 

Due to the symmetry of the problem we only evolve 
one octant of the grid. The spacetime is evolved using 
the NOK formulation and the excision methods described 
in 1^. The apparent horizon, if it can be found, is 
located using the horizon finder described in The 
gauge conditions used are "l-flog" slicing for the lapse 
and the gamma-driver condition given in equation (45) 
of 0. The event horizon is located using the code de- 
scribed in |33 |. 

In Fig. |31 we show the mass of the event and appar- 
ent horizons of the black hole found by measuring their 

1/2 

surface areas and using M = {A/16Try' and the ADM 
mass computed by the inital data solver for comparison. 
At early times the expected results are seen: the mass of 
the apparent horizon is less than the mass of the event 
horizon which is less than the total ADM mass of the 
spacetime. At late times large violations of the Hamil- 
tonian constraint develop, errors in the determination of 
the horizon become large and lead eventually to the crash 
of the code at around t = 0.455ms. Without excision the 
code crashes shortly after formation of the apparent hori- 
zon. 

We believe that the origin of these instabilities lies in 
the spacetime excision technique, and not the hydrody- 
namic excision described here. This can be seen much 
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FIG. 2: A simple shock propagating parallel to the x axis hits an excision boundary normal to the x axis (dotted line). As 
expected, the shock is completely absorbed. The solid line is the exact solution, the crosses are the numerical solution for the 
PPM scheme. Plots of the other reconstruction methods show qualitatively the same behaviour. The time shown is t — 0.5. 





FIG. 3: As in Fig.|21 a shock is propagated parallel to the x axis through an excision boundary. In this case a non-conservative 
excision boundary condition is used. This excision boundary condition, which is based on setting all variables in the first 
excised cell to very small ("atmosphere") values gives inaccurate but stable results. The circles represent the stable and 
accurate numerical result found using the excision boundary method described in this paper. 



more clearly in another test, where we used a slowly ro- 
tating neutron star (model Dl of 10]) and excise the 
spacetime variables and the matter variables at differ- 
ent times and compare the evolution of the 2-norm of 
the Hamiltonian constraint in the two cases. Note that 
it is not possible to first excise the spacetime variables 
and later the hydrodynamics because the hydrodynamics 
code is very sensitive to errors in the spacetime. There- 
fore we always first excise the hydrodynamical variables 
and then switch on the excision of the spacetime at a later 
time. We excise the hydrodynamic variables shortly af- 
ter an apparent horizon is found (which in this case is 
at a coordinate time of 0.546ms) and vary the coordinate 
time, toxj when the spacetime is first excised. 

Figure 1^1 shows the Hamiltonian constraint violation 
for simulations where the excision time tex is varied. For 
comparison, a simulation where neither the hydrodynam- 
ical variables nor the spacetime variables were excised, 
and a simulation where only the hydrodynamical vari- 
ables are excised, are also shown. 

Where the hydrodynamical or spacetime variables are 
excised the Hamiltonian constraint violation inside the 
excised region is not taken into account. At the point 
where either set of variables is excised the constraints 
are clearly meaningless within the excision region. We 



expect this to have no effect on the exterior spacetime as 
demonstrated numerically in, e.g., 0. The constraint vi- 
olations in the exterior spacetime clearly are meaningful, 
and we will be most interested in the behaviour shortly 
after either set of variables is excised 

It is clear that shortly after the time when the space- 
time variables are excised there is a sudden exponential 
increase in the violation of the constraints. This does not 
happen at the time when the excision of the hydrody- 
namical variables first takes place. We believe this indi- 
cates that the cause of the instabilities encountered is the 
method used to excise the spacetime variables. Whilst 
we cannot rule out instabilities from the hydrodynami- 
cal excision, the simulation where only the hydrodynam- 
ical variables are excised seems to lose accuracy due to 
classical slice stretching effects and not due to the high 
frequency oscillations which are instead observed when 
excising the spacetime variables. 



VI. CONCLUSION 

Long term evolutions involving black holes and matter 
will require excision boundary conditions for the matter 
fields. For hydrodynamics, where the fluid will gener- 
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FIG. 4: Michel solution compared with the numerical evolution for the ENO method. All other methods look very similar in 
the physical quantities. The graphs show the values for the x component of the velocity Vx, the density p, and the pressure p of 
the fluid. Each quantity is evolved until time t — 40 (in units of the background black hole mass), and the numerical values are 
compared with the exact solution (the solid line in the graphs). The x axis is the radial distance expressed in black hole units. 
The excision region is delimited by the vertical, dotted line. The last graph shows the convergence factor of all reconstruction 
methods, where the expected second-order convergence is seen for all methods except for PPM. 



ically form shock waves, the results of Lax and Wen- 
droff 16] and Hou and Lefloch 17] imply that a conser- 
vative total-variation stable scheme must be used, which 
presently means a HRSC scheme. Although errors in 
the hydrodynamics at the excision boundary should not 
propagate outside of the horizon and so should have no 
effect on the physics, such inconsistencies may lead to 
numerical instabilities. 

In this paper we have shown a simple method of provid- 
ing a consistent excision boundary condition for systems 
of equations in conservation law form evolved with HRSC 
schemes. Although we have only applied this technique 
to the Valencia formulation of relativistic hydrodynam- 
ics 0, Q, ^3 we expect the results to apply to any 
system where the conservation law form can be applied. 
The method is simple and works on a broad range of 
HRSC schemes. 

The tests of section [3 show that the method works 
for shocks, steady state cases, and dynamical evolutions. 
The shock tests are particularly simple and enlightening, 
showing that an incorrect boundary condition may lead 
to inconsistency or numerical instabilities. The Michel 



solution test shows that our boundary conditions are also 
suited to long term steady state evolutions with strong 
gradients. 

The dynamical collapse of a spherically symmetric 
TOV star shows the advantages and current limitations 
of our methods. The collapse to a black hole is accurately 
followed long past the formation of the horizon, the mass 
of which is computed to high accuracy. However, the 
simulation does not continue indefinitely, eventually fail- 
ing due to a numerical instability. Simulations including 
rotation show the same features and, by varying the ini- 
tial time at which excision is applied, there is evidence 
that the cause of the instabilities is the excision method 
applied to the spacetime variables, not the hydrodynam- 
ical method that is the focus of this paper. However 
we want to emphasise that although we believe that the 
simple lego excision boundary condition is the cause of 
the instabilities, the excision of the spacetime variables 
is nevertheless improving the length of the simlulations 
in the cases described here. 

It has been suggested that HRSC methods may im- 
prove spacetime evolutions for certain hyperbolic formu- 
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FIG. 5: Horizon masses for the collapse of a spherically symmetric neutron star. At early times the evolution is stable and 
accurate. The neutron star collapses smoothly into the horizon and the matter is excised. The horizon masses as measured by 
apparent and event horizon finders are very close to the total ADM mass of the initial spacetime, as expected. At late times 
instabilities lead to inaccuracies, as seen in the inset, before the simulation fails. 
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FIG. 6: The Hamiltonian constraint violation over time, plotted for different starting times tox of the spacetime excision. In 
every case (except the reference case where no excision is used) the hydrodynamical variables are excised as soon as the apparent 
horizon is found at t — 0.547ms. Instabilities at the excision boundary are believed to be the cause of the eventual failure of 
the simulation. Here it is clear that when the excision boundary condition is applied to the spacetime variables there is an 
immediate exponential growth in the constraint violation. In contrast, no growth is seen at the time when the hydrodynamical 
variables are excised. The instabilities which are indicated by the exponential growth of the constraint violations are clearly 
triggered by the spacetime excision. 
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lations of Einsteins Equations (e.g., [s^ l39l |40|'). If 

such a method were used, excision boundary conditions 
such as those presented in this paper should be easy to 
adapt. There are also other approaches to applying exci- 
sion boundary conditions that show great promise (e.g., 
^3)E3jE3|)- It however, not clear if these methods 
will work in situations where fundamental variables may 
become discontinuous. In such cases the methods pre- 
sented here could then be used for a matter field where 
the evolution equations could be written in conservation 
law form, and any alternative method could be used for 
the spacetime. 

At no point here we have addressed the question of 
the well-posedness of the system used, or any analytical 
proof of the stability of our method. The well-posedness 
of the simple Euler equations with arbitrary initial data 
in more than one dimension in the absence of gravity 
is not yet established; the most relevant current result 
( ^44] ) considers the one dimensional Euler equations on 
a plane-symmetric Gowdy spacetime. Given the current 
state of knowledge, and the successes in extending one- 
dimensional schemes that are known to be stable to more 
than one dimension, we expect that methods such as the 



one given here will be stable but are unable to give a 
proof at any level. 
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